home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_3.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.9 KB  |  81 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.3 (Taylor's Method of Order 4).
  9. % Section    9.4,    Taylor Series Method, Page 448
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements Taylor's method
  13.  
  14. % for solving the initial value problem
  15.  
  16. %     y' = f(t,y)    with    y(a) = y0
  17.  
  18. %    Define and store the function d1 = f(t,y)    
  19. % and formulas for the next three derivatives
  20. % of  f(t,y)  in the M-file  df.m 
  21. % function z = df(t,y)
  22. % d1 = (t - y)/2;
  23. % d2 = (2 - t + y)/4;
  24. % d3 = (- 2 + t - y)/8;
  25. % d4 = (2 - t + y)/16;
  26. % z = [d1 d2 d3 d4];
  27.  
  28. delete df.m
  29. diary df.m; disp('function z = df(t,y)');...
  30.             disp('d1 = (t - y)/2;');...
  31.             disp('d2 = (2 - t + y)/4;');...
  32.             disp('d3 = (- 2 + t - y)/8;');...
  33.             disp('d4 = (2 - t + y)/16;');...
  34.             disp('z = [d1 d2 d3 d4];');...
  35. diary off;
  36.  
  37. df(0,0); % Remark. df.m and taylor.m are used for Algorithm 9.3
  38. pause       % Press any key to continue.
  39.  
  40. clc;
  41.  
  42. %    Place the endpoints of [a,b] in  a  and  b.
  43.  
  44. %    Place the initial value y(a) in  ya.
  45.  
  46. %    Place the number of subintervals in  m.
  47.  
  48. a  = 0;
  49. b  = 3;
  50. ya = 1;
  51. m  = 12;
  52.  
  53. [T,Y] = taylor('df',a,b,ya,m);
  54. points = [T;Y];
  55.  
  56. pause    % Press any key to see the list of solution points.
  57.  
  58. clc;, clg;
  59. c = 0;
  60. d = 1.75;
  61. axis([a b c d]);...
  62. plot(T,Y,'g');...
  63. hold on;...
  64. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  65. if m<=30,
  66.   plot(T,Y,'or');
  67. end;...
  68. xlabel('t');...
  69. ylabel('y');...
  70. title('Taylor`s solution to y` = f(t,y)');...
  71. grid;...
  72. axis;...
  73. hold off;...
  74. shg; pause    % Press any key to continue.
  75.  
  76. Mx1 = 'Taylor`s solution to y` = f(t,y).';
  77. Mx2 = '     t(k)               y(k)';
  78. clc,echo off,diary output,...
  79. disp(''),disp(Mx1),...
  80. disp(''),disp(Mx2),disp(points'),diary off,echo on
  81.